home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / RUNGKUT.LBR / RUNKUT.FOR < prev    next >
Text File  |  1986-04-11  |  5KB  |  160 lines
  1. $NOFLOATCALLS
  2. $NODEBUG
  3. $STORAGE:2
  4. c*****************************************************************************
  5. c    NOTE: Some of the variables are not the same as in the article
  6. c          For example hadj = v, est=e(i+1), hfact=s, hlim = vlim.
  7. c*****************************************************************************
  8. c
  9. c
  10.     subroutine runkut(xa,ya,xb,neqn,tola,tolr,hstart,w,
  11.      &              imeth,ierror,icom,func)
  12. c
  13. c****    This is simply a front-end subroutine to split up the work
  14. c****    array specified by the user into smalller arrays
  15. c
  16.     implicit double precision (a-h,k,o-z)
  17.     external func
  18.     dimension icom(4),w(*)
  19.     common /epsil/eps
  20.  
  21.     call solver(xa,ya,xb,neqn,tola,tolr,hstart,ierror,imeth,
  22.      &            icom(1),icom(2),icom(3),icom(4),w(1),w(1+neqn),
  23.      &            w(1+neqn*2),w(1+neqn*3),w(1+neqn*4),func)
  24.     return
  25.     end
  26. c*****************************************************************************
  27. c
  28. c
  29.     subroutine solver(xa,ya,xb,neqn,tola,tolr,hstart,ierror,imeth,
  30.      &              ist,ichk,idef,iround,kt,est,yold,ynew,w,func)
  31.  
  32.     implicit double precision (a-h,k,o-z)
  33.     logical hflag
  34.     external func
  35.     dimension ya(neqn),w(13,neqn),kt(neqn)
  36.     dimension est(neqn),yold(neqn),ynew(neqn)
  37.     common /cons/hmin,hmaxl,hfact,hlim,power
  38.     common /epsil/eps
  39.  
  40. c****    If TOLA is set to zero, then a relative error test
  41. c****    If TOLR is set to zero, then an absolute error test
  42. c****    If neither are set to zero, then a mixed error test
  43. c****    If this is the first call to RUNKUT calculate the machine
  44. c****    epsilon and work out the constants depending on method used.
  45. c****    To ensure that the results will always be "safe", the value
  46. c****    EPS used by RUNKUT is actually 20 times the machine epsilon.
  47.  
  48.     if(ist.EQ.0) then
  49.       call caleps
  50.       call const(ist,imeth,idef)
  51.       eps=eps*20.d0
  52.     end if
  53.  
  54. c****    Check if the sum of TOLA and TOLR is greater than 10*EPS
  55. c****    If not, set TOLR to one million times EPS.
  56.  
  57.     ierror=0
  58.     if(ichk.GT.0)then
  59.       if(tola+tolr.LT.10.d0*eps)then
  60.         ierror=1
  61.         tolr=1.d06*eps
  62.       end if
  63.     end if
  64.     iround=0
  65.  
  66.     isig=dint((xb-xa)/dabs(xb-xa))
  67.     hold=hstart
  68.     x=xa
  69.     do 10 i=1,neqn
  70.        yold(i)=ya(i)
  71. 10    continue
  72.  
  73. c****    Set HMAX to the maximum value HMAXL set in CONST unless this
  74. c****    is greater than the interval width, in which case HMAX is
  75. c****    equal to the interval width.
  76.  
  77.     hmax=dmin1(hmaxl,dabs(xa-xb))
  78.  
  79. c****    call the runge-kutta solver using the current step size
  80. c****    to predict y(n+1)
  81.  
  82. 15    call rkp(x,ya,hold,neqn,imeth,kt,est,yold,ynew,w,func)
  83.  
  84. c****    calculate step size adjustment factor HADJ
  85. c****    then calculate HNEW using the smallest value of HADJ
  86.  
  87.     hadj=9999.d0
  88.     do 20 i=1,neqn
  89.        absest=dabs(est(i))
  90.        if(absest.EQ.0)then
  91.          hadjl=hlim
  92.        else
  93.          if(yold(i).EQ.0)then
  94.            hadjl=((tola+tolr*tolr)/absest)**power
  95.          else
  96.            hadjl=((tola+tolr*dabs(yold(i)))/absest)**power
  97.          end if
  98.        end if
  99.        if(hadjl.LT.hadj) hadj=hadjl
  100. 20    continue
  101.  
  102. c****    adjust the step size to HNEW using the calculated value of HADJ
  103. c****    unless HADJ is greater than HLIM, in which case use HLIM
  104. c****    If HNEW is larger than HMAX, choose HMAX as the new step size
  105. c****    If HNEW is less than HOLD/10. keep HNEW as HOLD/10 to avoid
  106. c****    excessively large swings in the step size
  107.  
  108.     hold1=dabs(hold)
  109.     hnew = dmax1(hold1/10.,
  110.      &             dmin1(hfact*hold1*(dmin1(hadj,hlim)),hmax))
  111.  
  112. c****    Check if HOLD is large enough compared to EPS to avoid
  113. c****    severe round-off errors. If HNEW is less than HMIN, the
  114. c****    problem has got too stiff or is discontinous-exit saving
  115. c****    the last points calculated.  If the last step was sucessful
  116. c****    let YOLD=YNEW and calculate YNEW using new step size. If it
  117. c****    was unsuccessful, recalculate YNEW using reduced step size.
  118. c****    If the HOLD was a reduced step size, then restrict HNEW to HOLD
  119. c****    If XB will be reached or exceeded,  exit after calculating
  120. c****    Y at XB.
  121.  
  122.     if(dabs(x).GT.0) then
  123.       if((hnew/(dabs(x)*18.d0)).LE.eps) iround=1
  124.     end if
  125.     if(hnew.GE.hmin) then
  126.       hnew=isig*hnew
  127.       if(hadj.GE.1.d0) then
  128.         if(isig*(x+hold-xb).LT.0.d0) then
  129.           x=x+hold
  130.           if(.not.hflag) hnew=hold
  131.           hflag=.true.
  132.           do 30 i=1,neqn
  133.          yold(i)=ynew(i)
  134. 30          continue
  135.           hold=hnew
  136.           goto 15
  137.         else
  138.           hstart=hnew
  139.           hold=xb-x
  140.           call rkp(x,ya,hold,neqn,imeth,kt,est,yold,ynew,w,func)
  141.           do 50 i=1,neqn
  142.          ya(i)=ynew(i)
  143. 50          continue
  144.         end if
  145.       else
  146.         hflag=.false.
  147.         hold=hnew
  148.         goto 15
  149.       end if
  150.     else
  151.       hnew=isig*hnew
  152.       ierror=ierror+2
  153.       xb=x
  154.       do 60 i=1,neqn
  155.          ya(i)=yold(i)
  156.  60      continue
  157.     end if
  158.     return
  159.     end
  160.